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FOREWORD 


This document represents Volume I of the final report 
on NASA Huntsville Contract entitled "Spin Vector Control 
for a Spinning Space Station". The report is prepared in 
two volumes: 

Volume I - User's Manual 
Volume II - Analytical Manual 


ii 



TABLE OF CONTENTS 

Page 

Foreword ....... il 

Table of Contents * H* 

Abstract i'V 

I Introduction. . i 

II Description of Russell's Method 3 

III Program Structure 17 

IV Description of Special Features 20 

A. Restrictions and Assumptions 20 

B. Movable Mass and Pendulums 22 

C. Attitude Control by Reaction Jets. . ‘ . . . 23 

D. Angular Momentum Control Devices 32 

E* Integration Method 41 

V Discussion of Results ...... 44 

A. Nutation Damping Pendulums 44 

B. Spin Control and Attitude Hold 46 

C. Propulsion Control 54 

D. Attitude Control by Reaction Jets 57 

VI References 60 

Appendix A. A**l 

Appendix B-l 

Appendix C "1 


iii 



ASTRACT 


This report constitutes the documentation of 
a technical study for developing an efficient com- 
puter program ror a restricted, rotational dynamics 
problem Involving multiple bodies connected by 
movable joints. The configuration is particularly 
applicable to a space station, system with a rotor 
producing a large spin vector. The computer' program 
is an effective tool for study of control of the 
spin vector and effects of the spin vector on other 


functions . 



I. INTRODUCTION 


The object of this study was to develop a means of investigating 
the dynamic behavior of a dual spin space station. In particular, 
the extremely large spin vector of such a system presents control . 
problems as well as having complex and often intuitively unexpected 
influence on the, response of "the system to internal and external 
forces. Emphasis was given to construction of an efficient computer 
program for simulating this system and functions relating to it. 

Basically the program is an application of Russell's momentum 
method of formulating and solving the equations of motion of a system 
of rigid bodies connected by movable joints. In the interest of 
efficiency the method is not applied in full generality. Rather, - 
the problem is restricted to configurations and constraints that 
are closely related to the actual space station systems that are of 
interest. Through choices gf parameters and optional features in 
the program, a sufficiently wide range of situations can still be 
accommodated. 

Primary features of the program are: 

1. ' Rotor-stator configuration with single degiee of freedom 

bearing joint. 

2. Complete generality in inertial parameters ’of both rotor 
and stator. 

- 3. Movable mass on rotor. 
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4. Nutation damping pendulums on rotor. 

5. CMGs on stator. 

6. Reaction jets on stator and rotor. 

7. Attitude control by reaction jets or CMGs. 

8. Transverse rate (nutation) removal by reaction jets, 

9. Spin control by reaction jets and/or torque motor. 

10. Flexibility of control functions by incorporation of 
control laws in subroutines. 

11. Gravity gradient effects. 

12. Docking capabilities. 

Although the program was exercised extensively in order to 
assure proper functioning, no exhaustive study of any problem 
or aspect of spin vector behavior was made. In the course of 
program check-out some interesting results were obtained and these 
are described. 
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IX. DESCRIPTION OF RUSSELL'S METHOD 


Equations of motion, pertaining to this problem are derived in 
detail in Appendix B. In this section the general development of these 

equations into a complete method of solution will be explained. 

The application of Newton's law to a configuration of rigid 

bodies' as given in Figure 1 results in equations of motion as 

follows . 

Derivative of Total Angular Momentum 


H 


+ T, 


EXT 


EXT 


2.1 


Derivative of Angular 'Momentum of Body 1 


E - + *0*1 + * 3 ~1 + *4*1 + 3 13 11 ? 3*1 

EXT 

2.2 

+ 3 14 * *4*1 + V* f 2*l 


Derivative of Angular Momentum of Pendulums 

- l 3 * f 13 t E 4 - V 4 - I 4 X F 14 


Since Body 3 and Body 4 are point masses, h 3 = h^ - 0 and 

• • 

= 0, for the pendulums then 




l 


* F 13 ! 


" 1*4 


K * F 14 


2.3 
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In the Euler method angular momenta are expressed in terms 
of angular rates. The derivatives of angular momenta then 
involve angular accelerations. The integration is carried out 
on the angular accelerations to obtain angular rates. At each 
step the angular momenta can be obtained from values of angular 
velocity and a knowledge of geometric and inertial properties. 

In Russell's momentum approach the integration is carried 
out directly on angular momenta. The angular rates at each 
step are obtained from the angular momenta. 

The momentum method was selected for this program. Rela- 
tively simple expressions for derivatives based on Equations 
2,1, 2.2, and 2.3 are used in the integration. The problem 

of obtaining values for angular velocities from values of 
angular momentum obtained by integration will be resolved later. 

In order to carry out the integration it is necessary that 
terms on the right hand side of the derivative expressions 
(Equations 2.1, 2.2 and 2.3) be determined entirely by infor- 

mation available at each point. That is, in order to proceed 
to the next time point it is necessary that the derivatives can 
be calculated. 

The terms T. and T_ would have contributions from 
1 EXT EXT 

several sources; 

1. Externally applied torques as given forcing functions 
ihich are known functions of time. 
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2. Externally applied control torques depending by some 
control law on geometry and/or rates. A current 
complete solution will be supplying geometric and 
rate information from which these torques can be 
obtained at each point. 

3. Externally applied torques depending on orientation 
(gravity gradient). Again, the current geometric 
information will permit calculation of these forces 
and torques at each point. 

The torques acting around the various hinge lines (Tq_^, 

T, T, ,) are all a result of control and/or friction 
1«**3 l-*-A 

functions. These torques are determined at each point through 
specific control laws from angle and angular rate data. 

The remaining terms (d 13 x F^.^, ^ x ^3 x 

etc.) in Equations 2.1, 2.2 and 2„3 are all torques due to 
interbody forces. These terras will be dealt with in. three 


steps : 


1. Forces will be expressed In terms. of center of mass 


accelerations (r^, r 9 , r^, r^) through Newton’s 


law. For example. 


r i-^3 = m 3 (R + 
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The cluster center of mass acceleration is found from 
the sum of all external forces 


R 


I & 


EXT 


then 

*13 = “3^3 * r 3 * K X ^ F EXT 

Through a change of variables the terms involving 
these accelerations will be replaced by terms 
involving velocities only. Continuing the same 
example, define 

h 3 “ x m 3* 3 

then 

* r mu 

H “ ^3 X “3*3 + T 13 " iT ^3 X ^ ? EXT 

The replacement of mass center accelerations by 
mass center* velocities introduces time derivatives 
of other terms such as h* d i3» etc ‘ 
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3. Finally, it is observed that velocities can be expressed 
in terms of geometric configuration and angular rates. 
Hence, the torques due to interbody forces are deter- 
mined by known external forces and torques and the 
angles and angle rates. 

In this modified form (see Appendix B, equations B_5„3, B5.5) the 
derivatives of angular momentum are determined by the known 
external and interbody torques, the known external forces and 
by angles and angle rates. A procedure to obtain current values 
for angular rates from current values of angular momentum is 
clearly essential. It Is also apparent that a concurrent inte- 
gration of angle rates to obtain geometric reatures such as 
gimbal angles and orientation is also necessary in the solution 
of this problem. The only remaining requirement for determina- 
tion of derivatives at each step is the means of deriving angle 
rates from angular momenta. 

This procedure will now be described. Angular momentum 
depends on geometry, inertial properties and angle rates. 

For the system of rigid bodies in Figure 1 the angular momentum 
expressions are: 

Total Angular Momentum 


H “ V 5 o + h‘\ + *o x Vo + *1 x Vi + ? 2 x V2 


+ r 3 x m 3 r 3 + r 4 x m^ 
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"Primed" Angular Momentum of Body 1 


h’ » I 1 *S 1 + m 2 ^ 2 x ^2 + m 3^13 + x r 3 + m 4^14 + ^ X 

"Primed" Angular Momentum of Body 3 and Body 4 

^ - h * ” 3 ^3 J % ~ K * Vi 

I 

The center of mass distance vectors, r Q , r^, . . . are 

expressible in terms of various vector parameters 3's, i’s 
(Figure 1 and appendix B equation B2,2). Consequently, the 
center of mass velocities can be expressed in terms of deriva- 
tives of the geometric parameters. (3 Q ^, 3 13 , ^ 3 , etc.). These 
derivatives are in turn each determined by the angular rate of 
the body in which the parameter is defined as a constant. For 
example, -3^ « x d Q1 . When all these substitutions have 
been made, the angular momenta will appear as sum£ of terms 
•with the following properties: 

1. All terms are of first degree in angular velocity. 

2. Aside from the angular velocity factors, no other 
factors in any term are explicitly time derivatives 
That is, the other factors can all be evaluated 
entirely from a knowledge of the inertial parameter^ 
and the current values of all angles. 
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It is necessary now to recognize that -the actual computation, 
must be 'done on scalar equations for the components of the 
vector equations that have been used up to this point. Scalar 
equations will be obtained for the following six components of 


angular 

momentum: 


1. 

H 

X 

x component of total H 


2. 

H 

-y 

y component of total H 


3. 

H 

z 

z component of total H 


4. 

h ix 

component of h| around- 

ulic uuvons trained x— axis 



of Body 1 


5. 

®3 

component of h^ around 

the unconstrained hinge 



line of Body 3 


6. 

g i 

component of h^ around 

the unconstrained hinge 



line of Body 4 


These six 

components of angular 

momentum are. expressed 


as functions of the following six components of angular velocity: 
1. angular velocity of Body 0 around the Body 0 


2 . 


to. 


Oy 


x-axis 

angular velocity of Body 0 around the Body 0 
y-axis 


3. . angular velocity of Body 0 around the Body 0 

z-axis 

4. angular velocity of Body 1 with respect to Body 0 

5 . fl angular velocity of Body 3 with respect to Body 1 

3 


11 



6, angular velocity of Body 4 with respect to Body 1 

The six scalar equations for the components of angular 
momentum in terms of the six angular velocity components are 
tediously long. However, since w T s occurred only :.n the first 
degree in the vector equations, these scalar equations are linear 
in the angle rate components. The six linear algebraic equations 
will be stated as: 



or h - Mw 

Then to is found from « = M ^h. The elements of M are 
calculated (in EMCALC subroutine) entirely from inertial pro- 
perties and a knowledge of the geometry of the system (angles 
and components of d^^, ^]_ 3 5 ^ 2> 3 The -angles and h are 

supplied by integration. Hence, the angular velocities can be 
obtained and the integration loop is closed. 
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The process can be summarized as follows: 


1 . 


At time t given h . u> , 0- , 9„ , 0, , calculate 

n n’ n In’ 3n’ 4n 


n 

2. Up date to t , } by integration obtain h - , 0. , 

Ht 1 n-r J- X 

Q 3 ’ °4 

n+1 n+1 

3. Set up M matrix from constants and 0 , 0„ , 

n+1 n+1 


0 


4 n+l’ 


,-l 


4, Solve a) = M *i h 

n+1 n+1 n+1 


5. Return to step 1 for t 


.Initially and the geometric configuration is given. 
Starting values for H^, and h'^ are computed from geo- 
metry and w in XDOT, the derivative calculation subroutine, 
which h* s to be called to supply the initial derivatives. 
Initial values for g^.and g^ are obtained by calculating M 
and using the initial values of to in the Last two equations 
of h = Mw. 

At each time point the orientation of Body 0 with respect 
to an inertial frame is updated by integration of quaternion 
parameters. The complete procedure is outlined in the flow 
chart of Figure 3 and diagrammed in Figure 4 in more detail. 
Figure 4 also indicates the effects of the moving mass or 
elevator (Body 2) and the CMGs, The moving mass contributes 
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to h and elements of the M matrix. The CMGs modify certain 

elements of M and modify H , H and H before the solution 

a y z 

for c*>. 
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FIGURE 3 PROGRAM FLOW CHART 
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Figure 4 Program Flow Diagram 
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III. PROGRAM STRUCTURE 

A general view of the main program will be presented with 
reference to the flow chart Figure 5 , 

The RCON section initiates the read-in of data for each run 
and stops the program when the last run has been completed. 

The INPUT section reads in and prints out all input data -for 
each run in accordance with the options that are included. 

The Initialization section resets initial values of certain 
variables, evaluates constants in accordance with input data and 
calculates initial values for integrated variables. 

In th • Integration section the angular momenta and angle 
variables -re updated in time by appropriate substitutions and use 
of the FOM > integration subroutine. Docking impulse contributions 
to angular momentum are also made at this point. 

The Quaternion section is logically a part of integration, 
but involves considerably more substitutions in the employment of 
quaternion parameters. The ROMS subroutine is used for integration. 
The result is an updating of the transformation from Body 0 to the 
inertial frame. 

In the CGNTM section, control functions by CMGs and reaction 
jets are activated in accordance with instruction flags in the input 
data. 

In the INVERT section the M matrix is constructed by calling 
the EMCALC subroutine. Contributions to the angular momentum by 
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terms due to the movable mass and CMGs are mede at this point. Then 
the system of equations h = Mw is solved for the angular rates by 
making necessary substitutions and calling the subroutine SYEQNS. 
There are two alternative routes here depending on whether the 
pendulums are used. If pendulums are included, a set of six alge- 
braic equations Is solved. Otherwise, a set of only four equations 
is solved. 

After the angular rates have been found, time is updated and 
the derivatives of the angular momentum variables are calculated 
by calling XDOT. The calculation of these derivatives completes 
the evaluation of all quantities pertaining to a point in time, and 
the results are printed out in the OUTPUT section. Again, the op- 
tions included will determine the printout structure. 

Before returning to the integration section for the next step, 
checks are made to determine if the parametric modifications required 
by docking must be made and also whether the run is finished. 
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I ENTER 



Figure 5 program Flow Chart 
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IV. DESCRIPTION OP SPECIAL FEATURES 


A. Restrictions and Assumptions 

The various features and options of the program are subject 
to a number of restrictions and assumptions which will be listed 
in this section. Besides providing a single convenient access 
to such information, this section will also give an overall view 
of the capabilities of the program. 

1. Body. Configurations 

Maximum of five bodies - the motion of one body 
(Body 2, the movable mass) is completely constrained. 

The program can be run withe ut the pendulums and/or 
without the movable mass. 

If either the rotor or the stator is to be eliminated., 
the mass and moments of inertia of the eliminated body 
may have to be approximated by small values. If these 
parameters are set to zero, the M matrix may be singular 
or poorly conditioned. 

2. Gimbal Constraints 

The joint between Body 0 and Body 1 is a single 
degree of freedom bearing. The bearing axis is 
parallel to the x-axis of Body 0 and coincident 
with the x-axis of Body 1. 

The hinges between Body 1 and the pendulums are 
single degree of freedom. The hinge lines must be 
parallel to the y-z plane, 

3. Translational Motion 

The orbit is circular with a radius specified in 
the input data. The orbital position is obtained 
by linearly Incrementing the orbit angle with time in 
accordance with a constant angular velocity derived 
from the radius. Variations in gravitational forces 
due to non-uniform or non~spherical earth are not 
accommodated. The orientation of the orbital plane 
and position with reference to the earth are not 
specified. 
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4. Attitude Control by Reaction Jet 

With the reaction jet attitude control scheme 
employed in this program the maximum angle of 
spin vector reorientation that may be commanded 
in a single maneuver is 60°. 

5. CMC 

A maximum of six CMGs can be located on Body 0. 

Each CMG may be zero DOF (reaction wheel) , 1 DOF 
or 2 DOF. 

6. Movable Mass (Body 2, elevator) 

Body 2 moves in a straight line in Body 1. The 
motion is entirely' constrained and is determined 
by the subroutines SCALC and SDCALC. 

7. Reaction Jets 

Control is limited to the x, y and z axes of Body 0 
and the x and y axes of Body 1. Torque on each 
axis is produced by two oppositely directed pure 
couples identical in position and strength. One 
jet control gain (for all four jets) and one dis- 
tance between jets (for both couples) can be speci- 
fied for each axis. 

8. Second Rotor 

The effects of a second rotor can be approximated 
by a reaction wheel on Body 0. The reaction wheel 
angular momentum would be adjusted to cancel the 
angular momentum of the first rotor Body 1. This 
approximation is in error due to lack of freedom 
in positioning and due to incorrect distribution of 
masses . 

9. Integration Step Size 

The time increment may be selected. The choice 
depends on rates encountered in the problem and 
accuracy desired. 
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B, Movable Mass and pendulums 

The movable mass Body 2 and the pendulums Body 3 and Body 4 
are Incorporated Into the system at the very basic level of the 
equations of motion. That is, the actual system equations of 
motion' are structured to include these bodies. Consequently, these 
bodies have already been considered in the detailed discussion of 
the method and in the derivation of the equations of motion in 
Appendix B. 

Terms due to Body 2, Body 3 and Body 4 appear profusely 
throughout the expressions for components of angular momentum 
derivatives and in the formation of elements of the M matrix. 

When these bodies are absent all of the terms due to them are 
rendered ineffective by setting in values of zero for appropriate 
mass and length parameters rather than being eliminated by pro- 
gramming . 
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C. Attitude Control by Reaction Jet 

The basic approach to the problem of attitude control by use 
of reaction jets will be introduced with reference to Figure 6 . 
Initially the spin vector (x-axis of space station) is aligned with 
x^., the x-axie of the inertial frame. It is desired to reorient 
the spin vector to a direction given by the vector c which here 
lies in the x^. - y plane at an angle of a from the x^. axis. The 
desired reorientation is accomplished by a simple rotation about 
the z axis. However, when a moderate torque is applied briefly to 
the z axis, the path taken by the spin vector will be the nutation 
cone A which misses the desired orientation by a wide margin. 

In principal a "brute force" technique could be employed to 

rotate the spin vector through the desired angle. A sufficiently 

s 

large torque applied about the z^ axis would ensure that the pre- ' 
dominant motion is essentially a rotation about the z ^ axis. As 
the direction c^ is approached, a large reverse torque is applied 
to remove the transverse rate. Transverse rates remaining after 
this cancellation can be removed by the normal propulsion control. 
Remaining direction error can be reduced by repeating this procedure. 

Instead of using this inefficient method the approach employed 
here minimizes fuel consumption by utilizing the nutational coning 
motion of the spin vector as part of the control policy. The 
desired orientation is approached by allowing the spin vector 
to traverse a segment of free body nutational motion (such as B 
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in Figure 6) resulting from application, of certain torques. At 
the termination of this path segment the spin vector direction 
is substantially closer to c . Residual transverse rates are 
removed by the normal reaction jet control. The process is 
repeated until the specified direction is approached as closely 
as required. 

The attitude control subroutine will now be discussed in 
more detail. Terms are defined in Table X.. 

A reorientation procedure will not be initiated unless 
transverse angular rates are below some prescribed minimum 
( |W0(2)J and |W0(3)| less than .0002). Transverse rate removal 
is given priority because these rates would interfere with the 
attitude correction process. 

If transverse rates are within tolerance, the spin vector 
direction is compared to the specified attitude direction vector 
CA. In this program the attitude test is made on attitude vector 
direction cosines • transformed to BODYO coordinates (CB) . If this 
difference is within tolerance ( <£L°), no attitude correction by 
reaction jets will be made. 

This correction process is open loop. The vector CB is 
not updated and no attitude error comparisons are made during the 
process. The sensed attitude error is used only to establish 
initial values and process parameters . 
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Process initialization consists ot storing tne current •J.'J.flO 
as TC, calculating normalized components GA1 and CA2 of the pro- 
jection of the attitude vector CA in the Y-Z plane of BODYO, and, 
calculating the torque magnitude AK required to produce the correction 
path. 

The required torque AK is estimated from simple approximate 
geometric and dynamic relationships. In Figure 7 the desired re- 
orientation of the spin vector from to x Q2 is a rotation through 
the angle ft . A nutation path of the H vector could pass near the 

directions x ni and x if the spin and transverse angular momenta 
U JL \J£ 

are related by sin 0/2 « . g^ rail - s - v . e - - s - e - . H trangver3e produced 

spin - 

by a torque AK acting for approximately 6 seconds, that is H transverse 
6 AK. The required torque is then approximated by 
AK * 1/6 Wl(l) * B0DY1I(1,1) ' sin 0/2. 

The discussion will be continued with reference to Figure 8 
a representation of the Y-Z plane of the TCI frame - a frame fixed 
in Inertial space coincident with the BODYO frame at the time that 
the attitude correction process is initiated. The projection of 
the desired direction vector on this plane is GAC. CA2 and CA3 
are direction cosines of CAC. The projection on the Y^^. - Z^,^ 
plane of the path which the spin vector direction takes due to appli- 
cation of torque AK for the tiime At^, is the arc B. The direction 
of torque AK is obtained from GA1 and CA2. 
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Z-TCI 



Figure 8 


Application of Application of 

gtarting torque Midpoint ending torque 



Figure 9 



The projected path B shows that the spin vector does pass 
close to the desired direction CAC, but it is necessary to apply 
& torque to stop the nutational motion near CAC* The estimate 
for the time of application of this second torque is obtained by 
sensing the midpoint of the arc. For this purpose define an angle 
0 between the TCI y-axis and the projection CAC. Note that cosy# 
and sincere already available as CA2 and CA3. In a frame obtained 
by rotating TCI through an angle (3 around the TCI x-axls the z 
component of the spin vector reaches its maximum excursion ZTC 
at the midpoint of the arc B. The time at which ZTC is sensed 
is the basis for timing the terminating torque near CAC. Note 
that the direction of the terminating torque is the same as for 
the initiating torque since the direction of travel of the spin 
vector projection has reversed in the course of traversing the 
arc B. 

At the time 2TMB-TMA when the terminating torque application 
is complete, the normal propulsion control is restored. After 
transverse angular rates are again removed, the process is re- 
peated to bring the spin vector successively closer to the 
specified direction. 

The reaction jets for removing transverse angular rates can 
be located on either the stationary or rotating body as provided 
for in PCON. The torques required for the correction process 
itself, however, must be obtained from jets located on the 
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Figure 9 Attitude Control by Reaction Jet Flow Chart 
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Table I Terminology in Attitude Control Routine 


AK. - Calculated estimate of torque for attitude correction 

AOJ (12) - Reaction jet lever arms 

CA2 ) - Direction cosines of projection of attitude 

CA3 J - vector in YZ plane of Body 0 at start of process 

CB(3) - Attitude vector in Body 0 frame 

FAT (8) - Force equivalents for attitude correction torques 

at reaction jet lever arms 
ICFA - Attitude correction process flag 

ICFB - Process initialization flag 

ICFC - Midpoint sense flag 

ICFD - Torque application flag 

TC - TIBO at start of process 

TMA\ 

TMB i - Times associated with application of torques 
TMC J 1 ' 

ZTC - Maximum excursion of. spin vector from "diameter" 
of correction path 

ZTCL - Previous value of ZTC} used in sensing extreme value 
of ZTC 
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stationary body. This restriction is necessary because jets on 
the rotating body will not in general be in the proper position 
at the time that torque is required. 

The attitude control by reaction jet feature was devised 
to provide large changes in direction of the spin vector. This 
restriction is consistent with the coarse tolerances, and crude 
approximations applied. It is assumed that fine adjustments 
to attitude and continuous accurate attitude control would be 
supplied by CMGs. 

The simple system simulated in this program is also limited 
in the reorientation angle that can be accommodated in a single 
command (a change in the vector CA). In order to avoid algebraic 
sign reversals that would disrupt the control process, the commanded 
change in spin vector direction must be less than 90°. A margin 
must also be provided for nutational path inaccuracies due to the 
crude estimates used to establish control process parameters. 
Accordingly, a maximum spin vector reorientation angle of 60° 
has been imposed for a command. Of course, attitude changes by 
larger angles can be accomplished by a succession of several 
commands . 

The simple reaction jet attitude control system presented 
here is intended to be an example of one of the problems that 
can be handled by the program. It is recognized that considerable 
Improvement could be made In this- system and that other systems 
may have superior performance. 
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D. Angular Momentum Control Devices 

The math model in this Section gives the angular momentum of 
the k UI CMG about the center of mass (c.m.) of the CMG, expressed 

in the frame. Control commands to the model will be gimbal 

• • 

rate commands 0 and 0, depending on the degrees of freedom 
possessed by the CMG. These rates are integrated to keep track of 
0 and 0. 

Total mass of the CMG is to be considered as part of body 0 

and located at the c.m. of the component. This point mass must be 

incorporated into the mass and inertias of body 0. The total angular 

momentum vector H of the CMG is to be expressed in the ' o _ frame and 

th 

included with all other such quantities on the o body to constitute 
the angular momentum used in Rnssell f s equations. 

Figure 10 shows the basic configuration necessary to describe 
the reaction wheel, one degree-of-freedom (DOF) CMG, or two DOF 
CMG. Considering the two DOF case, the gimbals are defined by 0 


and 9. 
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When 0 and 9 are equal to zero, the axes of the wheel are aligned 
to the zero gimbal frame (C k ). The frame Is then related to 
the body frame by the transformation £c 3 or by t * ie inverse 

relationship [ B o c 0 kJ ♦ 

Each momentum device will contribute a momentum vector consist- 
ing of two parts given as 


ok 


= f 


ok 


E , oj , 
ok o 


D-l 


The first part, £ q ^, is due to the commanded gimbal rates while the 
second part, E is due to the angular velocity of the body. 

As explained in section (II), the w's of the five main bodies are 
realized by solving the equation. 

w = M - h 


E , contributes to the M's and f . contributes to the h ? s. 
ok , ok 

Before attempting to derive the equations several variables 
need to be defined : 


1. 

B 

o 

Body frame for Body 0, 


2 . 

°ok 

Null gimbal frame for device 

in body 0. 

3. 


] Transformation from body 0 frame 
frame. 

to k null gimbal 
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4 . 


5 . 


6 . 


7 . 


8 . 

9 . 

LO. 


11 . 

12 . 

13 . 


14 . 


15 . 


16 . 

17 '. 


CVoJ In X erse of ^ G ok B oJ* 


ok 


tlx 

Total momentum vector of k device expressed 


. th _ 

In o frame. 


H,. 


H, 


H 


H 


th 


O) 


W.. 


Total momentum vector of k device expressed 
in null glmbal frame. 

Angular momentum of the inner gimbal Including 
mass of the wheel (C Qk frame). 

Angular momentum of the outer gimbal (C^) f^ame. 

Angular momentum of the wheel (C ) frame. 

ok 

Angular rate of the null gimbal frame. 


L w zj 

ok 

Prime Frame Principle axes of the outer gimbal. 
Double Prime Principle axes of the inner gimbal. 
fOG k C ok ] Transformation from null gimbal frame to 
outer gimbal frame. 

[ C ok 0G oJ Inverse ° f Kk C okl- 

[lG' ok C ok ] Transformation from null to inner gimbal 
frame. 

C G ok IG oJ Inverse o£ t IG ok C oJ* 

til 

H . Momentum of the k wheel. (Assume constant 
ok 

for CMGs ) . 
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18 . 


Inertia matrix of outer gimbal. 


19 # x Inertia matrix of inner gimbal including inertia 


of wheel. 


20. 

1 

0 

0 


Ka J - 

0 

C0 

S0 



L° 

-S0 

C0 

k 

21. 

C9 

S0S9 


-S6C0 

Kk c <J 

0 

C0 


S0 


S0 

•• 

-S0C0 


C0C0 

Degree-of -Freedom 

CMG 





For the 2 DOF CMG, the total momentum will consist of three 


parts 


H. ° H + H + H D-2 

k I o 

Each component of Equation D-2 will be derived separately. They 
then will be summed and the sum split into two parts; one part 
a function of the other part a function of the gimbal rates 
and wheel momentum. 

The angular velocity of the outer gimbal frane is related 
to the angular velocity of the body 0 by 


Mok m Kk c oJ t G ok B o3 W- 
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Therefore, the inertial rate of the outer gimbal is equal to the 
angular rate of the outer gimbal plus the gimbal rate 0 Q j c - 
The angular momentum of the outer frame is then given as 

E H o-Ut “ ^o-^ok { + ^ok} 

expressed in the prime frame and as 


w* = w og oJ ca* 

expressed in the null gimbal frame. When expanded, 

Wok - f°ok 0G ok] Wok { Kk°ok ] [c ok B o ] W + «ok}- D - 4 

The ttj. component of Equation D-2is obtained by applying a similar 

approach. Angular rate of the inner gimbal is given in terms of the 

body ra :e o> as 
o 

Mok - [ IG ok°ok ] [c ok B o ] [ "„ 3 - 
The inertial rate of the inner gimbal is given as 
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and 



When expanded and expressed in the null gimbal frame, H becomes 


[“ilk - t C ok 10 J Cl I 3 ok { Kk C ok 3 Kk’o 3 W 



D-5 


The third component H, of Equation D-2, can be obtained by- 
expressing the momentum of the wheel in the double primed system 
and then transforming it to the null gimbal frame. Let the 
momentum of the wheel be 
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since the spin vector of the wheel is always along the Z" axis. 
Therefore, 



This can be written as 



when expressed in the null gimbal frame. 

We now have the three components of expressed in the null 
gimbal frame. At this point we want to divide into the component 
due to and the part due to gimbal rates as given in Equation D-2, 

H, = E ,w + f . . 
ok ok o ok 
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therefore, when the three components of Equation D-2 are added to 

fch 

expressed in the o frame, the two components are 


E ok - tVoJ { [ C ok 0G oJ t;] ok [00 ok C ok ] [O ok B o ] 

D-7 

+ [ c ok iG ok] i>;P Kk c „J [Vo] } > 

F = [B C ] I [c OG ] [l ' ] 
ok L "o ok l ok ok o 




1 — 

•TS. 

1 


se 

+ C C ok IG oJ 

& G ok 0G ok3 

9 

0 

+ H 
z 

-ces0 

cec0 



One Degree-of-Freedom CMG 

Equation D-2 reduces to for the one degree-of- 

rreedom CMG, and the angle 0 no longer exists. By setting 0 and 
0 equal to zero. Equation D-5 reduces to 


t- H I-hk " 


+ 


Kk 0G <J 


{ Kk c oJ 


} 


ok 


[c , B ] [w ] 
ok o o 


D-9 
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and Equation D~6 reduces to 


M* - 


H 


z 


r i 


se 

o 



D-10 


Summing terms and expressing in the o 


th 


frame 




[Vot 1 { t°ok IG oJ hi’] [X0 ok c ok ] [o ok B 0 ]} 

D-1I 


- - 


r i 


0 


S0 

l?o C ck] { C C ok IG ok] hi'] Kk 0(! ok] 

0 

+ H 
z 

0 


0 


C0 


D-12 


Reaction Wheel 

For the reaction wheel, there are no terms that are a function 

A,1_ 

of the o body rates. Therefore, E^ is equal to zero. Also, since 

• 

0 and 0 are equal to zero, Equation D-12 reduces, to 

f = [B C ] H 
ok u o ok J z 
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E. Integration Method 

Components of angular moments, girabal angles and quaternion 
variables must be integrated in this program. In the Interest 
of obtaining computational efficiency some attention was directed 
toward selection of an appropriate routine for updating all 
these integrated variables. 

Three methods were considered. 

1. Predictor-Corrector 

2. Runge-Kutta 

\ 

3 . Adams 

A predictor-corrector method had been used in a previous 
program (Ref. 3 ). Some machine-time analysis showed that 

iterations of the lengthy calculations involved were very time 
consuming. Consequently, a non-iterative method was sought 
that would be accurate enough and yet avoid repetitions of 
the long calculations for each step. 

The Runge-Kutta method is non-iterative but does require 
calculation of derivatives at points between t n and t^ +1 in 
order to obtain the value of the variable at t In many 

applications the calculation of derivatives at the intermediate 
points is not difficult and the improved accuracy easily offsets 
this minor additional computation. 

In the rotational dynamic problems treated by this program, 
however, the calculation of derivatives is quite involved. The 
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derivatives of angular momentum depend on angular rates. These 
rates in turn are found by constructing and inverting the matrix 
which relates angular momenta to angular rates. In terms of this 
program a derivative calculation would require entering the 
INVERT BLOCK which calls the EMCALC and SYEQNS subroutines and 
finally calling the XDOT subroutine. Apparently most of the 
calculation for each step would have to be repeated for each 
intermediate derivative. It is doubtful then that the Runge- 
Kutta method can produce much improvement in speed over the 
iterative predictor-corrector method. 

The Adams method uses derivatives at t Q , t ... to 
integrate to t n+ ^ (Ref. 4 ). Since derivatives at these 
times have already been calculated for previous steps, it is 
necessary only to store them for future calculations. 

Comparisons were made between a 4th order Runge-Kutta and 
a 2nd order Adams method in a rotational dynamics problem. 

With the same value of At the Runge-Kutta was much more accurate 
due to the higher order of approximation. However, when At 
in the Adams program was reduced so that the accuracy of the 
two programs became comparable, the Adams program still showed 
a distinct advantage in speed. 

Some higher order Adams Integrations were also investi- 
gated. An increase in the order improves the accuracy without 
any appreciable increase in machine time. Serious disadvantages 
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arise though when the use of remote past derivatives produces a 
delay that leads to poor transient response or even instabilities 
in some situations. Many problems could be computed more 
effectively by a higher order Adams integration. However, since 
this program is to be applied to a wide variety of problems, it 
is advisable to remain with the simple second order Adams method. 
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V 


DISCUSSION OF RESULTS 


A. Nutation Damping Pendulums 

Preliminary nutation damping pendulum simulations were made 

4 

using moderately large pendulums disposed as follows on the rotor. 
Length - 10 ft. 

Hass - 100 slug 

Center of rotor to hinge distance - 100 ft. 

Even long runs (up to 1000 sec real time) with these large pendulums 
failed to show a perceptible decrease in transverse angular rates. 

Observation of small changes in transverse rates was hindered 
by the unsymmetric rotor. Angular rate variations introduced by this 
lack of symmetry tend to obscure small changes due to other effects . 

The possibility that stator pendulums would be more effective than 
rotor pendulums was considered. With appropriate parameter values and 
initial conditions this type of run was made using the same pendulums. 
Again no change in transverse rates was apparent. 

Reduction of At did not indicate that these results were due 
entirely to Instabilities introduced by computational errors. 

In order to exaggerate effects to observable levels, pendulums 
of these extreme proportions were simulated. 

Length - 200 ft. 

Maas - 3000 slug 

Rotor center to hinge distance - 300 ft. 

Under these conditions the pendulums on the rotor produced definite 
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Table 2 


Parameter Values for Computer Runs 


BOMASS - 
BODYOI =■ 

BIMASS = 

BODYlI = 

D01 « 0 
SP = .4 


4 

.6 x 10 

0 
0 

4.2 x 10 8 J 

.8 x 10 3 


1.4 x 10 
0 
0 


3.2 x 10 
0 


8 


6 x 10 7 0 

0 6 x 10 7 

0 0 


0 

0 

6 x 10 6 
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B. :5pin Control and Attitude Hold 

There are many techniques available that can be used to 
maintain the relative rotational rates between bodies O' and 
1 equal to a constant. The spin control loop that is con- 
sidered in this report is basically a one actuator control 
loop where the one actuator consists of a direct electric 
motor drive between the despun hub and the rotor. A control 

law which drives the torque, motor and adaptively compensates 

\ 

for frictional torques has been programmed in Subroutine 
TORKO.1 and is given by the following equations : 

T KOTOR ^ S 2 + ^ 

£ 3 So ) 

DESIRED 


5-1 

5-2- 


with the initial conditions 


£(o> - - g 2 n 1 


5-3 


It is easy to show that the equilibrium occurs when 

** ant * that this solution is independent of any 
frictional torques, A typical sequence of events follows. 
Assume that the frictional torque suddenly increases. The 
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result Is a decrease in the spin rate and a corresponding 
increase in £, which is coupled via equation 5-1 to give an 
increase in T M0T0R resulting in an increase in to the 
desired value. In this manner changes in friction are adaptively 
compensated. 

A simple approximate analysis can be used to show how the 
gains g£ and can be determined. Assuming that bodies 0 and 
1 are sufficiently decoupled so that we can write Euler's 
equation as : 


« « 

hAx " < t motor " £rictl<ra) 

X 0x e 0x = " ^MOTOR " frictio11 ) 
recalling that 


lx 


fi l + 


Ox 


5-6 


Substitution of equations (5-4) and (5-5) into 9^. =» + 9^ 

yields 


h 

i 


(T. 


MOTOR 


- friction) (T 


MOTOR 


- friction) 


Ox 


'lx 


5-7 
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Differentiating (5-7) and assuming 


0 we get 



’’’MOTOR T MOTOR 

I T 

Ox lx 


5-8 


Differentiating equation (5-1) and substituting into (5-8) 
gives with some rearranging 


“l + 


I 0x + I lx 


l3ox' I 'lx 


g 2 fi l 


I_ + I, 
Ox lx 


LVix 


“ g 3% 


’’Ox ^ I lx q 

I_ I, g 3 DESIRED 
Ox lx -> 


5-9 


which is the' desired result. 

Equation (5-9) is a simple linear second order differential 
equation. It is now an easy matter to choose the gains and 
gg. to give any desired dynamic response. Figure 11 shows a 
typical time response of The desired value of was .4 rad/ 

sec while the actual value was .39 rad/sec when the problem was 
started. After an initial transient the actual and desired 
are approximately equal. 

Once the relative angular rates oetween coaxes x ana u 
are stabilized to some desired value there remains the problem 
of orienting body 0 in inertial space. A simple CMC control 
law was used for this task. Two single degree of freedom 
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SPIN ANGULAR VELOCITY OMEGA1 RAD/ SEC 



0 . 


1.000+01 2.000+01 3.000+01 4.000+01 5.00 


OEGA1 


VS 


TIME fN SfzC. 


Figure 11 


Spin Control Time Response 
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CMGs were positioned on body 0 so that there respective 
momentum vectors in the null position are opposite each 
other and perpendicular to the momentum vector of the 
body as is shown in the sketch below: 


^BODY* ~*TOT 


H. 


CMG2 


H 


GMGl 


A ^BODY 



H CMG2 ' H CMG1 


Clearly,, 


^OT = ^ODY + H CMG 


5-10 


H CMG = 2H CMGl 8in 9 CMG 


5-11 


Differentiation of equation (5-11) yields 


H CHG ” 2H CMGl COS 9 CMG 8 CMG 


5-12 


As a control law for equation 5-12 let 


9 CMG " S 4 6 0x : S 5 9 0 x 


5-13 
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Using the uncoupled form of Euler's equation and assuming that 
the only torque acting on body 0 is that due to t.he CMGs the 
following equation results 

C 0x°0x = 2H CMG COS 9 CMG ^ g 4°0x + ^Ox' 1 5-14 

For si.iall CMG angles cos9 ~ 1 and equation (5-14) degenerates 

CMC? 

to 


1 0x 9 0x “ 2H CMG s 4 9 0x + S 5°0x 5-15 

Equat Lon (5-15) can be used to select the constants and 
so that the dynamic response of 0^ is satisfactory. Using 
value.-; for g^, g g , I Qx is shown below in Table 3 

a com) uter run was initiated whose results are shown in figure 12 


Table 3 


Constant 

Value 

8 4 

-1 

g 5 

-1 

h cmgi 

7.0 x 10^ ft-lb-sec 

v 

1.4 x L0 7 slug ft 2 


Figure 12 shows how the CMGs repositions the angle 0^ from 
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T!B0(3 ) '2) 


1,000+01 2.000+01 3.000+01 4. 000+01 5.00 

vs time: in sec . 


Figure 12 Attitude Hold Time Response 
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and initial offset value of 9 a ,175 radians to the desired 

ux 

value of 0 Qx = 0 radians. This is simply a representative 
run showing how CMGs can be used to reposition or reorient 
the attitude of body 0. 
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G. Propulsion Control 

A demonstration o£ transverse angular rate reduction by use of 
reaction jets was made fox the standard space station configuration 
(Table 2 ). Control gains were - 2 x 10 5 and jet couple arm lengths 
were 200 ft. 

With jets on the stator an initial transverse rate magnitude 
was reduced by a factor of four in approximately 30 sec. (Figure 13 ) 
The transverse rate reduction by stator. jets was in good agreement 
with results predicted by the formula given in'tha User’s Manual 
description of PCON. 

The effect of an identical set of rotor reaction jets is shown 
In Figure 14 , The transverse rate "magnitude was reduced by a factor 
of four in 100 sec. The rotor jet action is slower because the rotor 
has jets on only one transverse axis, which due to rotation is at 
times not oriented properly for effective control. In terms of 
total impulse required for a given reduction in transverse rates, 
however, the effectiveness of rotor and stator jets was about the 
same. 

In Figure 13 and Figure 14 it is seen that the stator acquired 
a small spin axis rate (W0(1)) during the control action, A correspond- 
ing change in spin rate of the rotor also occurred. These rates 
could have bean corrected by spin axis jets. Spin, axis control was 
not applied in these runs -in order that transverse control comparisons 
would not be influenced by spin axis control 'effects. 
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RAD/ SEC 



o. 2.000+01 a. 000+01 e.ooo+oi 8.000+01 i.oo 


TIME IN SEC 


DELTAT = .5 SEC 

TIBOI(l,l) - TIBOI(2,2) - TIBOI(3,3) * 1.0 
WO(l) « O WO (2) = .02 RAD/SEC WO(3) = 0 
OMEGA1 » .4 RAD/SEC (INITIAL) 


AOJ.(l) = 10 
AOJ (2) a 200 
AOJ (3) = 200 
A1J(1) - 10 
AlJ (2) = 10 


CGAINO(I) =0 
CGAINO (2) = - 10T 
CGAIN0(3) = - 10 5 
CGAINl(l) = 0 
CGAIN1(2) = 0 


SEE ALSO 
TABLE 2 


(INITIAL) 


2 


Figure 13 


PROPULSION CONTRA- 
JETS ON BODY 0 


RAD/ SEC 


2.500- 02 HJjjjjTFJ 

fete 

2.000-02 

1.500- 02 
1 . 000-02 



- 2 . 000-02 


-2.500-02 


0 . 


2.000+01 

4.000+01 6.000+01 

8.000+01 

1.00 


TIME IN SEC 




DELIA! - .5 SEC 

TIBOI(l,I) - TIBOX(2,2) « TIBOI(3,3) - 1.0 
WO(l) = 0 WO(2) ■= .02 SAD/SEC WO(3) = 0 


AOJ(l) 

AO J (2 ) 

AOJ (3) 

A1J(1) 

A1JC2) = 200 FT CGAIN1('2') = - 10' 


i 

WO (2) ■= .02 RAD/SEC 

.4 RAD/SEC 

(INITIAL) 

10 


CGATNO(l) =■ O 

10 


CGAIN0(2) = O 

10 


CGAIN0(3) - 0 

10 


CGAINl(l) = 0 

200 

FT 

CGAINim = - 


SEE ALSO 
TABLE 2 


(INITIAL) 


Figure 14 PROPULSION CONTROL 
JETS ON BODY L 
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D. Attitude Control by Reaction Jet's 

The performance of the reaction jet attitude control 
function is indicated in Figure 15 . Initially the spin .axis 
is aligned with the inertial x-axis (TIB0(1,1) = 1.0, 

TIB0(2,1) = 0, TIBO(3,l) -0.) The direction to which the 
spin axis ia to be reoriented has the directions in inertial 
space CA(1) = .5, CA(2) =■ - .612, CA(3) « .612, 

The first attitude correction segment was completed in 
about 54 sec. The stator x-axis was rotated to a direction 
given by the direction cosines TIB0(1,1) - .63, TIB0(2,1) = - .6, 
TIB0(3,1) = .5. From 54 sec to 90 sec residual transverse rates 
were removed, and the spin axis direction changed only slightly. 

From 90 sec to 144 sec the second correction process occurre 
The spin axis arrived at a direction in Inertial space given by 
the direction cosines TIB0(1,1) = .52, TIB0(2,1) = - .62, 
TIB0(3,1) = .6. Subsequent removal of transverse rates occurred 
until about 174 sec. 

A third correction maneuver was initiated, buc was not com- 
pleted by the end of the run. 

At the end of the 200 sec run the new direction of the spin 
axis was within 3° of the specified reorientation direction. 
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DIRECTION COSINES OF SPIN AXIS IN INERTIAL FRAME 



•0. 4. 000+01 8.000+01 1.200+02 


TIME IN SEC 

DELTAT = 1.0 SEC 

TIBOI (I, 1) » TIBOI (2,2) - TIBOI(3,3) = 1.0 
WO(M) = 0. (INITIAL) 

OMEGA1 = .4 RAD/SEC (INITIAL) 

CA(1) = .5 CA(2) = - .612 CA(3) » .612 


AOJ(l) = 10 
AO J (2 ) = 200 FT 
AOJ(3) = 200 FT 
A1J(1) = 10 
AlJ (2) - 10 


CGAINO(l) =0 5 

CGAIN0(2) = • 2 x 10^ 

CGAINO (3) = - 2 x 10 

CGAINl(l) = 0 

CGAINl(2) = 0 


SEE ALSO TABLE 2 


1.600+02 2.00 


Figure 15 REACTION JET ATTITUDE CONTRO] 
:>b 


Instability, An initial transverse angular velocity -aagnitude of 
,02 rad/sec increased to more than .03 rad/sec in 150 sec. Pendu- 
lums on the stator had a definite nutation damping effect, The 
transverse rate magnitude decreased from .02 rad/sec to ,018 rad/sac 
in 150 sec of real time. 

These studies lead to the following conclusions: 

1, Reasonably sized pendulum nutation dampers will not be 
very effective, 

2, This program is not an. efficient means of studying 
pendulum nutation dampers unless they are of im- 
practically large proportions-, 

3, The inertial configuration proposed is unstable with 
nutation damping pendulums on the rotor. 
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APPENDIX A 


THK QUATERNION SOLUTION OF THE BODY TO INERTIAL COORDINATE TRANSFORMATION 


Consider two sets o£ coordinates systems whoso axes are denoted 
by X., Y . , Z. and X, , Y, , Z . Initially the two sots of axes are 
coincident. ^By Euler ' s theorem any rotation of the X fa , Y^, system 
may be expressed as a rotation through some angle, about some axis. 
Suppose then that the X^, Y. , Z set is rotated through an angle (the 
positive direction of rotation Is given by the right hand rule) about 
an axis whose direction cosines are cosa, cos/3 and cos 7 relative to 
the X., Y i , Z, set. Let x^, y fe , z ^ be the X^, Y b , Z^ coordinates of 
a vector and let* x , y^, z^ be the X , Y , Z coordinates of the same 
vector. The transformation matrix A of Equation A-l is required. 


A 


- - 
X I 


*b 


*= A 

y b 

z. 


Z, 

I 


b 


It can be shown 
l-2s 2 f b 2 a 

2 (s 2 cac/3 

+ s f c f c7) 

2 u 

2(cac7s — 

- ■ f • I •» 


that A is given fay: 


2 (s 2 - cac/3 


-sf c|c7) 
l-2s 2 | s 2 /3 


2(s 2 | c/3c7 

+ S ~ c ^ ca) 


(A-l) 


2 (cac7s 2 ^ 

+ * f c f 

2(s 2 | c/?c7 

- 8 f C 2 


(A-2) 


A. C. Robinson, "On the Use of Quaternions in Simulation oi 
Rigid Body Motion", TR58-17, Wright Air Development Center. 
Wright -Patterson AFB, Ohio, March 1960, pp 5-7. 
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In order to introduce the quaternion make the following sub- 
stitutions : 

e Q = c f, e L = cctB e 2 - c/3s e 3 - c7s ^ (A-3) 

The A matrix becomes 


A = 


The quaternion was invented by Hamilton in 1843. The quater- 
nion q is composed of four parts. 


2.2 2 2 
8 0 + a l ' e 2 - e 3 


Z(e l e 2 ' e 0 e 3^ 


2(e l e 3 + V2> 


2 < e l e 2 + V9 


2 2 2 2 

e 0 " e l + e 2 - e 3 2(e 2 e 3 - Vl> - 


2 < e l e 3 - e 0 e 2 


2( e 2 e 3 + e o e i) 


2 2 2 2 
e 0 " 6 1 " e 2 + e 3 


q *• q 0 + 1 + J q 2 + k q 3 


where q , q^, q^ and q^ are real numbers and 1, j, and k are hyper- 
imaginary numbers satisfying the conditions 


2 2 2 
i D j = k « 


ij « - jl = k 
jk - - kj = I 
ki *= - ik = j 


The conjugate of the quaternion q is : 


q* * q - i q. - j q_ - k q . The magnitude of q is defined as 

kl - (q*q) 


(A-4) 


A-2 < 



Now that the quaternion is defined its use in coordinate 
transformations will be developed. Define a quaternion by 


v r ix b + jy b + k v 


Now consider the quaternion product q V^q* where 

q - e + ie + je + ke . e , e , e , e are defined by 

Equation A-3. Then q V^q* “ (e Q + ie^ + je 2 + ke 3 > (ix b + jy b +'kz b ) 

(e - ie. - je - ke ). When this expression is expanded 
using the defined quaternion operations the following is 
obtained 

q V B q* - [(e 2 + ~ 4 - 4) *b + 2(e i e 2 • e 0 e 3 ) y b + 2 ^l e 3 - Z b] 5 

+ [2( ei e 2 + e Q e 2 ) ^ + (e 2 - e 2 + e 2 - e 2 ) y b + 2(e 2 e 3 - e^) ajj (A-5) 


+ [ 2 ( 6^3 - e Q e 2 ) ^ + 2 ^ + y b + <e 2 - e 2 - e 2 + e 2 ) *Jk 


Now if the A matrix of Equation A-4 i 3 substituted into Equation' A-l anc 
the indicated multiplication is carried out, the following expression 
is obtained 


J 

K 


^0 + 6 1 " ®2 " e 3^ *b + 2 ( ' e l e 2 " *0 e 3^ y b + 2 < - e l 6 3 + e 0 ® 2 ^ % 

y i 


2 (e l e 2 + e 0 e 3 ) *b + (£ 0 " e l ■ t 4 “ e 3 } y b + 2<e 2 e 3 “ e 0 e l ) Z b 

Z i 


2 <e i e 3 " e 0 e 2 ) *b + 2<S 2 e 3 + e 0 e l ) y b + (£ 0 - 4 ~ 4 + e 3 ) Z b 


(A- 6 ) 


The x, , y , x of Equation A -6 are identical to the coefficients of i, 
j and k in Equation A-5. Hence, the quaternion multiplication qV^q* 
accomplishes the same transformation as Equation A-l as long as 
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e Q « cos - cosa sin “ cos/? sin “ cosT sin 

Define the quaternion by 
V I “ ix l + jy i + kz i* 


Then 


Vj - qV B q* - qV g q* 


Now consider two successive rotations of the X^, Y, y system. 
Let the first rotation be defined by q. . Then suppose the 
X., Y, , Z, system is again rotated. Let this rotation be defined . 
by q 2 and b x^, y£, z^ denote the final position of the moving axes. 
Note chat q. must be referenced to X^., 
enced to X fe , Y b , Z fe . Let x, , y. . z 


the Xj, Y. 
vector. 


T» Z T* X b* A b* 
Then define 


S’ 


Y' 

b* 


4* 


Yj, Zj. and q„ must be refer- 
’ y b> V *b*/b i 

coordinates of th< 


b be 
.e same 


V I “ ix i + Jy i + kz i 

\ - i! s + 3y b + k2 b 

V B ' 1X b + iV b + kZ b‘ 


So that 

v ] " ’i’A* 

V b“-^W 


or 


v i- 




* = 


q l q 2 V B (q l q 2 ) ** 


A~4 



Thus the quaternion which defines the transformation from the 
X^, Y^, system to the X^,, Y^, Z i system is q^. 


Suppose the angular velocity of the rotating coordinate 
frame is known in terras of the components of angular velocity 
along the rotating coordinates . Let P, Q, R be the X- , Y,, Z, 
components of the angular velocity. Let q(t) define the relaiic. 
ship between the rotating and the fixed coordinates at time t. 
During time t to t + At suppose the moving set undergoes a rota- 
tion of Am about an axis which has direction cosines cosa, cos/3, 
cosY with respect to the moving frame . The quaternion q(t + At) 
which defines the transformation between the fixed and moving 
frames at time t + At is required. 

Now the transformation between the rotating frame at time 
t and the rotating frame at' time t + At is 

q A => cos A^ + sin Ay (icosa + jcosjS + kcosT). 


Since q(t) and q^ define successive rotations. 


q(t + At) must be 


q(t + A) « q(t)q^. 


Then 


q(t + At) - q(t) = q(t)q A - q(t) 
or 

"q(t + At) - q(t) = q(t) [_cos ^ + sin ^ (icosa + jcos/3 + kcosT)] - q(t) 
As At -*-0, Am-^0 so cos Ay and sin Ay — 

q(t 4- At) - q(t) = q(t) + q(tj,~ • (icosa + jcos/3 + kcosY) - q(t) 
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or 


« lim q(t + A^) " -. I (*• - ) . o iim ^ ~p - (icos a + jcos/3 + kcos7). 

At-*-0 At-^-0 

q(t) - \ q(t) (iP + JQ + kR). 

So that 

e Q + 1 1 + i 2 + k 3 = \ (e 0 + le l * je 2 + ke J } (1P + JQ + ^ 

“ \ ie 0 P + je Q Q + ke Q R - e 1 P + ke^ 

- je^R - ke 2 P - e 2 Q + ie 2 R + “ ie^Q - e^R 

" 2 < ' e l P + ®2^ + e 3 R ^ + 2 ^0 P + 6 2 R ” e 3^' ,t 
+ | (e Q Q - e L R + e 3 P)j + (e Q R + e^ - e 2 P)k 

Then the differential equations for the four components of the quater- 
nion are: 

S 0 = ‘ I (e l P + e 2 Q + e 3 R) 

* J < e 0 P + e 2 R " “ e 3« 

®2 “ I (e 0 Q " S 1 R + e 3 P) 

| (QqR + - e 2 P) 
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0, 


Thus, only four linear differential equations need be solved 
to find the components of the quaternion. The A matrix can then 
be determined from the components of the quaternion. Referring to 
Equation A-4 


The initial conditions are a (o) = 1, e^o) = e (o) = e <o) = 
since the moving and inertial system are coincident at t = 0, 


*11 “ 
*12 - 
*13 = 
*21 = 
*22 “ 


' *32 “ 
*33 = 


2<e l e 2 - V3> 

2 ( e l e 3 + e o e; 2' t 
2^ e l e 2 + e 0 e 3^ 


e 


2 

3 


2 < e 2 e 3 - VP 
2(0^3 - 


2 ( e 2 e 3 + e o <S l^ 

2 2 2 , 

8 0 - e l - e 2 + 0 


2 

3 


The A matrix computed from the quaternion components will be 
orthogonal if and only if the quaternion has a magnitude of unity. • 
From the definition of the components of q(t) the -nagnitude of q(t) 
should always be unity. However, the numerical integration of 
the differential equations for q(t) introduce errors which cause 
•,the magnitude of q(t) to differ from unity. q(t) may be kept at 
unit magnitude by normalization at each intergration step. Nor- 
malization is accomplished by dividing each component of q(t) by 
the square root of the sum of the squares of the components. 



APPENDIX B 


DERIVATION OF DYNAMICS EQUATIONS 

1. The Space Station Configuration. 

The equations of motion governing the dynamics of a dual spin 
space station will he derived by the momentum techniques as exploltec 
by Russell ■> [l,2] for an N-body spacecraft. The space station con- 
figuration is shown on Figure 1 and consists of a despun (stator) 
body and a rotating (rotor) body which contains two inverted pen- 
dulums and an elevator. Active control will be employed between the 
rotor and the .stator to maintain the correct spin rate of the rotor 
and overcome the frictional torques between the rotor and the 

stator. The attitude and nutation control will be achieved by 

\ 

control moment gyros in the stator and a suitable combination 
of thrusters on both the rotor and the stator 

The variables used in the formulation of the dynamics of the 
space station are defined on Figure 1, and with a few exceptions 
is basically the symbolism adopted by Russell loc. cit. The majo: 
assumption involved in the formulation of the dynamics of the 

1 

space station is that the two pendulums and elevator are point 
masses, the inference being that they do not have any moments 
.of inertia^ In these cases the bas ic\ r otational dynamics for e ach 
point, mass body is replaced by a .static torque balance which 
results in considerable simplification of the dynamics. 


B-l 


2. Center of Mass Equations 

From the definition of the center of mass it follows that 


Vo + Vl + m 2 r 2 + »3 r 3 + ra 4 r 4 ' 0 


B2.1 


The geometry of the multiple connected rigid body (see Figure 1) 
yields : 


r l " r o + d oi 


X — IT *4™ d $ 

2 0 01 2 


r 3 ~ r 0 + ^01 + *13 + ^3 


B2.2 


*4 = *0 + *01 + *14 + h 


Using equations B2.2 in equation B2.1 determines r^ in terms of 
the multiple body vectors as follows. 


Vo + ra l<*0 + *01 } + m 2 c? 0 + *01 + + m 3 (? 0 + *01 + *13 + ^3> 

+ + V + 5 ia + -0 = °? 


Defining m = m Q + m^ + + m^ + the^ above result can be re- 


duced to: 



mr n + < m " m 0> d 01 


m 2^2 + m 3 


(d i3 +i 3 ) +m 4 (d u + lj = 0. B2.3 


Substituting this result into equations B2.2 yield* the local center 
of mass vectors from the total center of mass position in terms of 
the multiple body vectors. 


r 


1 



^ 2^2 ^3 ^ d 13 + ^ 3 ^ 

m m 


m 4 (d 14 + h 


®0 , ... m 2. 2 
r 2 " ~ d 01 + (1 ' - sr 


m 3 (d 13 + ^> 


m 4 (d 14 + V' 


m 


B2.4 


- "o a 01 m 2^2 n "3,,, m 4< 3 14 + ^> 

r 3 " — - — + (1 - — > Cd i3 + 1 3 > - r 


3 <5 X3^ 3 > + (1 .^ ( , it+i) 


4 m 


14 


3. Rectilinear? Dynamics 

The rectilinear dynamics for the five body cluster depicted 
on Figt re 1 are determined by the following vector differential 
equatic ns (Newton ' s Law) . 
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m 0 (R + ;„> “ 5 F 0 « + F l-^0 

m l< R + r ]_) “ ^ F 10 + F 0-^l + F 2->-l + F 3 -HI + F 4 •►I 

1^2 ( R + " F l-*-2 B3.1 

m 3 (R + r 3 ) - F^ 3 

m 4’( R + *4^ = F l<-4 

The external forces acting on the stator (0 body) and the rotor 
(1st body) are denoted by Fq^ and respectively. The inter- 
body forces are represented by F > y, which is interpreted as the 
the force due to the body acting on the £- th body. Since every 

force has an equal and opposite reaction this yields; the skew 

\ 

relationship. 


V ss - F 

fx->-y y — 


The vector sum of equations B3.1 yields 


** + <Vo + Vl + m 2*2 + m 3*3 + V4> “ ? F 0« + 


*10 
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which reduces to 


irR 


? F 0« + ^ F l/S 


B3.2 


by virtue of the center of mass definition (B2,l). 

The primary purpose of equations B3.1 is that they are used 
to eliminate the unknown interbody forces which appear in the 
rotational dynamics. 

4. Rotational Dynamics 

The rotational dynamics for the five body cluster are 
governed by the following vector differential equations. 


Sa x F. + d n1 
a Oa Oa 01 


F l~0 


+ T 


l-*-0 


B4.1 


fi 


1 


- ^h/3 * + 3 13 * F 3~l + 3 u x F 4*l 

+-? a x ? 2 -o-l + + T 3-»-l + 


B4.2 


The angular momentum vectors ror the 0 body and the 1 body 

are denoted by h Q and h 1 respectively. Using the same notation 
as for the Interbody forces, the interbody torques are represented' 
by with the equivalent relationship T^.^, = - 

The torques due to gravity gradient will be represented by 
appropriate external forces acting on the bodies of the cluster. 
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Since we are dealing with point masses for the inverted 
pendulums, the rotational dynamics degenerate into a static 
balance of torques . 


T a x F 
A l-*-3 3 l-»2 


T h ii x F 

a 1-*-4 4 l-*-4 


B4.3 


The total angular momentum for the cluster of five bodies 
is defined by 


H - £ 0 + fi l + J o * Vo + i i x Vl + *2 x Vz 

0 • 

+ ; 3 X m 3 f 3 + ; 4 X m 4 i 4 


and is governed by the differential (Euler's) 'equation 

h - E(i 0 1 i 0a ) x f 0a + tyh + V 51 hf> B4 - 5 

The rotational dynamics of the space station will- be 
governed by the unconstrained components of the vector equations 
B3.1.(3), B4.2, B4.3 and B4.5. Physically the dynamics of the 
stator and rotor are given by equations B4.2 and B4.5j and 
equation B4.5 is chosen rather than equation B4.1 since it 
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avoids some of the constraint problems between the stator and the 
rotor. The effects of the elevator on the dynamics of the space 
station are inherent in the interbody forces of equation B3.1(3). 
The unconstrained component of this equation can be used to 
control the motion of the elevation, however, for this study we 
shall ignore the elevator control problem and determine the effects 
of prescribed elevator motions on the stability of the space 
station. Finally, the effects of the inverted pendulums are 
determined by equation B4.3. 

For these equations the unknown interbody forces have to 
be eliminated and the equations for the unconstrained components 
recast into normal form suitable for numerical integration. This 
basically is the momentum technique. 

5 . The Momentum Method 

The unknown Interbody forces are eliminated from the appro- 
priate .rotational dynamic equations by means of the rectilinear 
dynamic equations B3.1. Treating equation B4.2 first, and 
eliminating the interbody torques and in addition 

to the interbody forces yields 



~h xi 2 + “3 (d 13 + X *3 “ m 4 (d 14 + h X r 4 

-v ‘ 3? 

" [ m 2^2 + m 3 (d 13 + ^3 } + m 4 (d 14 + X ( 5"ir 
+ x + T oi B5, 
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If we define 


h' “ h x + m 2 H 2 xi 2 + m 3 (d 13 + x ^ + ^4) x * 4 

which represents the angular momentum of the rotor, elevator and 
inverted pendulums about the hinge point between the rotor and 
stator, then equation B5.1 can be put into normal form as 

dhf n. - "ft • 1 j i 

dT " * *2 + B 3< d 13 + V * r 3 + m 4 Cd 14 “ V * r 4 


- [™2^2 + m 3^ d 13 + + m 4^ d 14 + ^4^ x 


0 a 

m 


v_M) + Yjt 


+ + ^ a l/3 x F 10 + T 0-l 


B5.3 


Equation B4.3(l) for the inverted pendulum reduces to 


M • A 


T l-^3 ml 3 * V* + r 3> 


i 3 x m 3 i 3 +i 3 K )(X)5 0a + £ V 


Defining 


S 3 “ ^3 X *3*3 


B5.4 
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seduces the pendulum equation to 


p-.i 

dt 3 


x m 3 r 3 



m 


x 


m 


p o« + 


+ 



B5.5 


Similarly by defining 


S 4 


l. X 

4 


m, r, 
4 4 


B5.6 


the second pendulum equation becomes 


dh! 
4 

dt 


m 


£>. x m. r, - & , x 
4 4 4 4 m 




a 




+ T 


1—4 


B5.7 


Therefore, the unconstrained components of equations B5.3, B5.5 

and B5.7 together with equation B4.5 represents the rotational 

dynamics for the space station. The total angular momentum 

fctl 

differential equation will be expressed in the stator (0 body) 
coordinate system, therefore equation B4.5 becomes 


dH 

dt 


+ cu x ' H 


£<; 


+ a n ) x F n + 


a 


j a i 
0 1 


) X F. 


If we substitute for r Q by equation 2.3 and r^ by equation 2.4 
then the above reduces to 
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dH 


tt V T (m_m 0 ^ 5 m 3 (3 13 + ^3 ) ®4^14 + 'V 


+ oi x H - Z-, L" 


dt 0 


m 01 m m 


m 


. - n = ,v r 5j ^2 V J 13 + ^3 > V J 14 + V 

+ a 0a J x F Qa + 4-/ [ m d 0 i ~ m - m - m 


+ a 10^ x F 10 


B5.8 


Similarly the angular momentum and the inverted pendulum equations 
will be expressed in the rotor (body 1) coordinate system* so that 
equations B5.3, B5.5 and B5.7 become 


dh 

dt 


i ^ 

“ +a) i x “ m 2 ^2 ** 2 + m 3 ( ^13 + X *3 + ffl 4 (5 14 + X *4 


- [”^2 + ” 3 < a i3 + i 3 ) + * 4 (a 14 + i 4 )] X <12— + T, 


+ £a lfl x r 1(J + T 0 ^ 


B5.9 


dh * 
3 


V 1 - 


dt* + *1 x £ 3 " *3 * Vs - l S * ST <5 F 0« + + T l-~3 B5 - W 


and 


31T + Sj X h| - i 4 x m 4 ; 4 ; i 4 X ^ (E F 0a + + T^ 4 B5.ll 
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To develop the unconstrained components of these equations the 
relative orientations and klnematrics for the cluster of bodies 
must be defined. 

Define j «= + “3^13 + ^ 3 ) + m 4^14 + ^u) 

Then B5.8 can be written as: 


dH 

dt 


+(D, 


x H ° 


1 

m 


m “ m 0 - - 

~ d A1 - j) X 2jFr, 

m 01 J/ a 0« 



(m 0 d 01 


- 5) x 


Ef. . 

0 1 P 


+ 


Ea rt x F„ + 
a Oar 0a 


Ea^x 


■10 


B5.12 


and B5.9 can be written as: 


dH| + x h^ m 2^2 x *2 + ra 3^13 + ^3^ x *3 + m 4 ( ^14 ' + K') X *1 


Jt 


“ ra ^ ( 5 P 0 « + ^ F l/P + ^ a l/3 X *10 + T 0-^l 


which can be written 


dh‘ + ^1 x h[ = ?i x [- m 2 l z - m 3 (d 13 + I 3 ) - m 4 (d 14 + l^] 
dt 


" m ^ ^ E f q or + Ef^) + Ea^ x + T Q ^ B5.13 
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Apply equations B5.10 and B5.ll to find the components along 
the unconstrained hinge lines of the pendulums. 

Unconstrained component along axis 


dt ^ 3 y S 3 y + h 2z B 3z^ + S 3y ^“lz h 3x “ ^IxNz^ + s 3z ^lx h 3y ** ^ly^x^ 
= m 3 s 3 y <i 3z i 3x " ^3x r 3 z) + “3 S 3z ^3x r 3y " Ky^x* 

+ 5 T a 3 y {^<* 0 * + V ' V" V iU9 l + F 0 z COs9 l + F l*>} 


, m 3 

+ rrT S 3z 


{4,<v oae i + E o* 8ln9 i + V - V F <* + r ix4 B5 - 14 


Unconstrained component along s^ axis 


ar < h 4y a 4y + h 4z S 4z ) + a 4y (w lz h 4x " W lx h 4^ + S 4z (w lx h 4y ' to ly h 4z ) 


~ m 4 s 4y ^4z r 4x “ K **4^ + m 4 S 4z ( ^4x r 4y “ ^4y r 4x ) 

+ r a 4jr {4z<V + F l*> - ^4x ( V 1 ” 9 ! + F 0z“ s6 l + F lz>} 

* IT s 4z'{ 4 X <V 0s9 1 + F 0z al " 6 l + V - V F 0* + F l* ) } B5 - 15 
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APPENDIX C 


GRAVITY GRADIENT EFFECTS 

In this Appendix a derivation will first be outlined of 
the expression for the gravity gradient moment about the mass 
center of a single vehicle segment. The total gravity gradient 
moment about the mass center of a cluster is regarded as the 
sum of the individual moments thus computed for the segments, 
plus a- moment due to differential forces acting at the mass 
centers of those segments. The differential force expression 
is then derived under the assumption that the entire mass of 
a segment is' located at its center of mass, and results from 
the difference in position within the gravitational field of 
the segment’s mass center and that of the cluster as a' whole i 
1. The Gravity .Gradient Moment on a Single Segment 

Consider Figure C-l, with R q being the radiu.3 vector from 
the center of the earth and r, that from the center of mass 
of the body to an increment of mass dm. 

The force on dm is then 


IF 


- u dm R 



(C-l) 


C-l 




Figure C-l 


The resulting moment is 


dM = r x dF 


- u dm rxR 



Making the substitutions R = R + r and 

o 


R 



3 . 





we can write the moment as 


C-2 



dM 


u dm r x (R q + r) 


1 + 


2R • r 
o 

R 3 


3 

2 


. (C-3) 


Since 


r x (R q + r) = r x R q 


2R * r 

i + - _ + £_ 

3 2 

R R 

o o 


2 

2 V 2 


and 


1 - 


3R * r 
o 


for r « R , 
o 


we can integrate (C-3) over the entire body to get the approxi- 
mate expression 


M 


/- 

Bj 


u r x R 


R 


3R * x 

° * 1 ^ !dm. (C-4) 


Expanding and noting that J" rdm - 0, 

Bj 


(C-4) can be reduced to the form 


M = - 


3 u 


R' 


J 


(R x r) (r - R ) dm. 


o Bj 


(C-3) 


If this vector equation is written in matrix form, referring 
all vectors to the body-fixed B^ -frame, and the indicated 
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integration carried ou:, we arrive at the desired expression 


(3x1) , (3x3) (3x3) (3x1) 

M .-V E I. [R ] (C-6) 

R 5 3 ° 

o 

where [r 1 is the mat: ix expression of R in the B. -frame, and 
R is the cross product matrix obtained therefrom: 



and Ij is the convent] >nal moment of inertia matrix of Body J 

in terms of the B -frame. Equation (C-6) is used in the gravity 

1 

gradient subroutine in the form shown, 

2. Differential Force on Body J 

Consider Figure C-2, where r^ denotes the vector from the 
mass center of the cluster to that of Body J, having mass m.. 

The total force on Body J is 



CLUSTER 



MASS 


Figure 02 

Following a procedure similar to that of paragraph 1 above, 
the approximation 


u m 


F j 58 


■i ) n 


R + x 
° J 


3R * r 


R 


1 R j 

° ) 


(07) 


is obtained, valid for R q . 


u >a,S. 

Since - — J* - ° is the gravitational force on B. that would 
TC J 


occur if r j were zero, the desired differential force is simply 


“"li- 3R o ■ r l - , 
5F . s-*Mr, — ” l R_}-. 


'J 


R' 


3 I r j 


(08) 


R 


The matrix equivalent of Equation (08) is used in the Gravity 


05 



Gradient Subroutine. 

The derivation of components of the vector R q in the frames 
of Body 0 and Body 1 will be described with reference to Figure 



Figure C-3 


The translational motion is restricted to a circular orbit 

in the X^-Y^. plane. The origin of the I frame coincides with 

the center of the earth. The orbit angle is measured from the 

+ Y axis. The vector R is easily expressed in the I frame as 
I o 



- l R ol 

I Frame 


-80 

o 

c0 

o 

0 
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Application of the inverse of the TIBO transformation will then 

give the components in the Body 0 frame. The components of a 

unit vector R /I R | are designated as DB(1), DB(2) and DB(3) 
o o 

at the beginning of the subroutine. For the Body 1 calculations 
these components are expressed in Body 1 coordinates through the 
A transformation. 
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